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Abstract 

Molecules with complex internal structure in time-dependent periodic potentials are studied 
by using short Rubinstein-Duke model polymers as an example. We extend our earlier work 
on transport in stochastically varying potentials to cover also deterministic potential switching 
mechanisms, energetic efficiency and non-uniform charge distributions. We also use currents in 
the non-equilibrium steady state to identify the dominating mechanisms that lead to polymer 
transportation and analyze the evolution of the macroscopic state (e.g., total and head-to- head 
lengths) of the polymers. Several numerical methods are used to solve the master equations and 
nonlinear optimization problems. The dominating transport mechanisms are found via graph 
optimization methods. The results show that small changes in the molecule structure and the 
environment variables can lead to large increases of the drift. The drift and the coherence can be 
amplified by using deterministic flashing potentials and customized polymer charge distributions. 
Identifying the dominating transport mechanism by graph analysis tools is found to give insight in 
how the molecule is transported by the ratchet effect. 
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I. INTRODUCTION 



Theoretical research on Brownian motors and the ratchet effect has flourished since the 
early 1990 's [l-4|. Most studies have been limited to simple systems with one or two coupled 



particles, whereas research of more complicated systems has escalated in recent years |5Hl5|. 
Due to the increased complexity of the models with internal structure, numerical methods 
play a more important part. This is due to the fact that the ratchet effect occurs in a far 
from equilibrium environment and only simple model systems can be analyzed exactly (see, 

e.g., m^). 

The internal structure is an important aspect for many real- life molecular motors (e.g., 
the well-studied kinesin il8|). For systems with non-homogeneous potentials, internal states 
usually play a more important part than in "traditional" transport driven by biased ex- 
ternal forces (such as a constant electric field). Even single particle systems based on the 
ratchet effect have been shown to display many phenomena, of which the current inversion 
phenomenon is one of the most interesting. Current inversions are found to be rather com- 
mon and can usually be generated by tuning of variables (e.g., diffusion constant, friction, 
potential shape and/or period) jl. [19I26I. In view of this it is reasonable to assume that 
systems with internal dynamics possess even more surprising properties, and knowledge of 
the correlation between internal states and transport would enable artificial engineering of 
the molecules and to boost wanted properties such as the velocity. An intriguing possibility 
considered in this Article is the control of electrophoresis j^^l by modifying the internal 
charge distribution of the molecule. Because of the large number of parameters and dif- 
ferent models, it is hard to compare results from different works and form any universal 
rules about the current or energetic properties for the ratchet effect. Things get even more 
complicated for complex molecules, for which the results are even more model dependent. 
Therefore we think that it is necessary to at least develop some general methodology for 
how to systematically study and monitor the behavior of these systems. This is indeed one 
of the key themes of this work. 

We have recently studied polymers using the Rubinstein-Duke (RD) model in time- 
dependent periodic potentials [28]. The RD model |29| is a good prototype of a complex 
molecule since the size of a linear polymer can be easily varied, it is strongly correlated, 
and the model has been actively studied for two decades 30-3^. There has also been in- 
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terest towards polymers as Brownian motors recently j^, Ig], [s^I • In Ref . [28[ we presented a 
general "toolbox" based on the numerical solutions for master equations and found current 
inversions for the RD model in the flashing ratchet and traveling potentials. In this Arti- 
cle we extend our work and methodology by considering the efficiency, different potential 
time-dependency schemes, non-homogeneous charge distributions and the dynamics of the 
internal states leading to the macroscopic transport properties. We formulate the opera- 
tors and master equations that are then solved with suitable numerical tools that fall into 
areas of linear algebra, integration, optimization and graph analysis. Due to the nature of 
the ratchet effect, most observables that we are interested in (such as drifts and conforma- 
tional changes) are very small. Therefore we find that a discrete space model that allows 
numerically exact solutions provides a very useful framework in this context. 

It is found that, like in many other studies on the ratchet effect before, varying certain 
model parameters has a large effect on the velocity, coherence and energetic efficiency. We 
take this aspect a step further by doing multiple parameter optimization for the RD model 
in order to maximize the steady state drift. If the internal states and the movement of 
the polymer are tightly correlated (such as in the RD model), changing the parameters 
increases the importance of some molecule conformations over the others. We demonstrate 
this by comparing the expected values for certain characteristic macroscopic properties for 
polymers, such as head-to-head and total length. We also identify and compare the most 
important microscopic conformations of the polymer that are responsible for the currents in 
different situations. 

This Article is organized as follows. In Section II a mathematical framework and notations 
are defined and in Section III we go through the numerical methodology. In Sections IV 
and V we present our results for transport properties and examine their microscopic origin. 
In the Appendix, the operator algebra involved is discussed in detail. Our conclusions are 
given in Section VI. 



II. MODEL 



We study the transport of the RD polymer 29| and its modification, the free motion 
(FM) polymer 28|], in temporally and spatially changing driving potentials. Essentially the 
RD model consists of connected Markovian random walkers {reptons) in continuous time 
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(see Fig. [T]). Each repton carries a charge that interacts with the potentiaL The model was 
originally developed to study the reptation process of the polymer in a restrictive medium 
(gel). However, in the context of this study, the model is primarily used as a good prototype 
of a molecule with a large number of internal states. To study the importance of the bulk 
motion, the assumption of the reptation can be relaxed, which results in the FM model. The 
complexity of the polymers can be increased by considering arbitrary charges of the reptons. 
In the following Section, formal definitions of the model are given for the implementation 
for numerical computations. Readers not interested in the technical details may skip this 
part and proceed to Section B. 



A. Stochastic generator and operators 



Consider a one-dimensional discrete Markovian random process in continuous time 



After the transition rates between all the allowed states in the system are given (el- 
ements Hi j), the stochastic matrix H can be defined. For molecular motors, this matrix 
includes all the internal conformations and spatial positions of the molecule in the potential 
3-411- In the case of the Markovian stochastically driven potential, it also includes the 
states of the external potential. We consider systems with stochastic {type 1) and deter- 
ministic potential switching schemes with sudden {type 2) and smooth {type 3) switching. 
The potential V{x,t) is assumed to be L and T periodic in space and time (for stochastic 
switching T is the expectation value). The stochastic matrix for the polymer dynamics is 



H 



.1=1 



y \ j=i / 



Us 



for the type 1 and 



N-2 



1=1 I y \ i=l / . 

for the deterministic case, where Tg is the expected lifetime of the potential Vg, and q G 
the repton charges. The switching of the potential is assumed to be cyclic, i.e., Vi — )■ V2 — >■ 
■ ■ ■ — 7- — )■ Vi. The operators rig and hg create transitions between the potential states, 
and the operators A, B and M determine the dynamics of the head {A and B) and middle 
(M) reptons. The detailed definitions of these operators are given in the Appendix (see 
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FIG. 1: (Color online) Illustration of a single configuration of the six repton (N = 6) polymer in 
the external potential V{x,t) with available moves shown by arrows. The moves with blue (gray) 
arrows are only available in the FM model. 



also Ref. 



28|). In Fig. [T] we have an illustration of the six repton polymer in one of its 



configurations. We fix the direction of the motion such that up arrows indicate the positive 
direction and vice versa. 



The type 2 operator now becomes 



Hi 



te[0,Ti) 

te[Ti,Ti+T2) 



t e 



and for the type 3 we choose V{x,t) = Knax(a;) sin^(7rt/T). The type 2 and 3 potentials 
are more reasonable for artificial molecular motors that have external driving mechanisms 
(e.g., electric potential), whereas the type 1 occurs most likely in nature (e.g., ATP driven 
motors). After the generator is defined, the dynamics is given by the master equation 
dP{t)/dt = H{t)P{t), where the elements of the probability vector P{t) include all the 
individual states y of the system. The stationary state Pgtat for the type 1 generator means 
that HPstat = and for types 2 and 3 that Pstat(^) = -Pstat(^ + T). From Pgtat, all expected 
values, such as the drift f , can be computed. The effective diffusion coefficient D^q is found 
by solving another stationary state for the diffusion equation (see Ref. |28|]). After solving 
V and -Deff, the Peclet number can be computed from Pe = 

Although the drift and Peclet number are the most studied properties, they tell nothing 
about the internal dynamics of the molecule. More specific measurements are needed. For 
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a discrete model, individual states, which we call microstates, can be bunched together to 
define the macrostates. Within the operator formalism, the general form of such macrostate 
operator is 

k yeFk 

where Uk is the corresponding value of the macrostate (e.g., the polymer length), riy a 
microstate operator, and Fk a (large) collection of microstates. For the RD-type model, 
there are 3^~^ microstates, for which the operators have the form 

N-l 

1=1 

where the function g{y,i) defines the state {A, B or 0) of each bond i between the reptons 
i and z + 1. We define the following four macrostate operators for the RD-type model: the 
zero-bond count (number of 0-bonds), the kink count (number of AB or BA bond pairs), the 
head-to-head length (distance between first and last repton) and the total length (maximum 
distance between two reptons). The head-to- head and total lengths are calculated in the 
potential direction (the only spatial direction for the one- dimensional model) and for the fully 
accumulated polymer they both are zero. The corresponding operator definitions of these 
observables are found in the Appendix. Separating the head-to-head and the total length 
is important since the polymer can take a U-shape. For example, for the configuration in 
Fig. [1] the values for Uk of these operators would be 2 for zero-bonds (formed by reptons 
1-3), 1 for kinks (reptons 3-5), 1 for the head-to- head length and 2 for the total length. 



B. Selection of the rates 



Despite the large number of studies with discrete state Brownian motors, the importance 
of choosing the rates Hij has not got much attention. By demanding the local detailed 
balance (no net currents in equilibrium), the usual choices for the rates are 42|, |43 | 



min {l, e(^J-^')/'=BT j (Metropohs) 



,{Ej-Ei)/2ki 



(exponential) 



^1 ^ g(i?,-E,)/fcsTj 1 (Kawasaki) 

where F sets the time-scale and l/ksT is the Boltzmann factor. Both of these constants 
and the lattice constant, are set to 1 in this paper. All three definitions lead to the required 



Pi = exp{—Ei)/Z distribution in equilibrium, but generate the different Icinds of dynamics 
when apphed to ratchet systems (far from equihbrium) such that the microstate energy Ei 
contains the potentiaL To demonstrate this, we have plotted in Fig. [2] the stationary state 
drift and diffusion coefficient of the 8-repton RD polymer in the type 1 flashing ratchet 
and traveling potential (model parameters are listed in the beginning of the section IV). 
Although all three curves for flashing and traveling potentials share a similar shape, the 
scales are different and large differences can be seen in the limit where the temporal period 
T — )■ 0. Being fast and simple, the Metropolis form is usually the favorite choice for the 
rates. But especially with ratchet systems it can be a poor choice since it does not take into 
account the slope of the downhill moves (rate being limited to 1) that is important for the 
dynamics. This is also true for the Kawasaki form, since it is basically just a smoothened 
Metropolis function. Since there is no single correct choice for the rates (based on theory), 
the selection must be made on experimental or model specific grounds. Only exponential 
(in flashing ratchet) and Metropolis (in traveling potential) dynamics lead to zero drift in 
this limit, which is a physically more realistic situation and is also consistent with the single 
Brownian particle model j^^. Therefore we choose these rates in this study. 



C. Non-uniform charge distributions 



The usual assumption in the studies concerning polymer transport is that all monomers 
are identical, i.e., they carry identical charge and mass. We relax this assumption and study 
the effect of the non-uniform charge distributions along reptons. Previous works on the 
RD model have considered some aspects of this. In Ref. 
(i.e., one charged head repton) was considered and in Ref. 
comes to the drift velocity all charge distributions are equivalent in small fields (i.e., linear 
response regime) 65|]. In Ref. |46| it was noticed that the drift in constant field depends 



a magnetophoresis model 



45| it was shown that when it 



strongly on the position of the charged repton within the polymer and in Ref. [35| non- 
homogeneously charged RD polymers in large fields were studied. Recent study of the dimer 
in the periodic potential show that if the connected particles are non-identical, directed drift 
can be generated even in the symmetric potential |7|. 

We want to find the best possible charge distributions q for the RD and FM polymers by 
finding the largest possible drifts. This leads to a mult i- dimensional, nonlinear constrained 
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FIG. 2: (Color online) The effect of the jump rate scheme. The drift velocity and the diffusion 
coefficient of the 8-repton RD polymer as a function of the temporal period T in the type 1 
flashing ratchet (left) and the traveling potential (right), with exponential (blue dash-dotted lines), 
Metropolis (black solid lines), and Kawasaki (red dashed lines) rates. 

optimization problem with constrains coming from the charges Qi. We choose Xli •?« ~ Q 
and Qi > 0, where the first constraint simply sets the total charge corresponding to an 
uniformly charged polymer and the second one fixes the sign of the charges. The optimal 
charge distribution gives some (indirect) information about the polymer conformations and 



V). Lastly we 
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49|. 



reptons that dominate the transport (i.e., have the largest impact on the dri: 
note that optimization has been carried out for some single particle systems 

III. NUMERICAL METHODS 



A. Network analysis 

The stochastic matrix H can be also treated as a graph with vertices (states) and edges 
(transitions) that can be analyzed to gain more detailed information of the transport process, 
as described in this section. Graphs and statistical physics have a long history due to the 
close similarities between stochastic systems and electric circuits, and in the seminal work 



of Schnakenberg in 70's many important results between these two were presented 



50j (see 
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5l| for some recent developments). Most of the works on this subject deal with the relations 
between steady state, rates, probability fluxes and entropy. We are however interested in 
finding the optimal paths within a current graph, which has not gained interest within 

previous works. Such ideas have however risen in other disciplines such as microbi olog y 



52| . In the following, only basic knowledge of the graph theory is expected (see, e.g., |53|). 
For simplicity, we consider only type 1 scheme where the time-dependency of the stationary 
state does not need to be explicitly dealt with thus making the numerical computations 
easier. 

After the stationary state Pstat of H is found, the net currents (edge weights) between 
the states can be computed. In addition to the stationary state and stochastic generator, 
we make use of the matrix H^^^^ which carries the information about the direction and 
magnitude for the transitions in H . For the RD-type model, the elements of if^'sn ^re ±a/N 
for all polymer state transitions, where a is the lattice constant and the factor 1/N results 
from the center-of-mass motion, and zero for the potential state transitions (present only 



for the type 1 system). For additional details, see Ref. j28|]. The graphs G and G'^^s" ^^j.^ 



then formed as follows. Let Wij = Psta.t{i)Hij — Pstat (j)-f^i,t V2,j. If Wij > 0, there exists 
a directed edge z — j in G and G^^^^ with weights Wij and w^'^^f^ = -^Hf^^Wij. With the 
sign in front of the weights w^lf^, one chooses the direction of interest of the transport (see 
below). The weights Wij are probability flows and the weights wfl" are mean displacement 
flows in the stationary state. 

Let 7fc be a path — )■ 22 — ^ ■ ■ ■ — ^ ^fc in the graph with 7^ iy'ix 7^ y (the path is 

non- intersecting) and 7fc(j) = ij- We then look for the path(s) 

fc-i ^ 

'^^'^yiT^ik{i+i)nk{i) ='■ ^^^fhk), (3) 

7fc ^—f K, 7fc 

1=1 

with X being Wij or w^^^^. The resulting path computed with Wi^j contains transitions that 
lead to the largest mean probability flow and we denote it by 7^. Similarly with one 
gets the path with the largest mean probability flow and we denote it by 7^^^". We call 
these paths the dominating processes. The function / is known as the target function. If the 
system is closed (periodic), the process must eventually return to its starting state and a 
cycle is formed, in which case 7fc = 71. Since the potentials we study are indeed periodic, we 
concentrate on closed systems from now on. For the cycle 7^^^", the target function defines 
the mean cycle velocity, i.e., Vc = Yl^Zi ^y^^^,yi^ where yi = 7fc^°(i). 



Whether there is a difference between 7 and depends on the details of the system. 
It may turn out that 7 only includes transitions that are not responsible for the directed 
molecule transportation, but instead results from the non-transporting diffusive motion. 
Formally this means that Yli=i = 0) which we call a stationary process, as the 

net transport for the cycle is zero. This is indeed typical for the ratchet transport, since the 
molecule spends most of its time near the minima of the potential, being unable to move 
until the suitable state of the potential and molecule conformation is reached. Therefore 7'^'sn 
carries more interesting information as it takes into account the directions and magnitude 
of the moves. If the path has a property Y2'i=i ^^^^{i+i) 'y^.{i) 0; we call it a transporting 
process. It is not guaranteed that the dominating process is a transporting process in either 
case. 

In the literature, the problem in Eq. ([3]) for cycles is known as the optimum cycle ratio 
problem (see, e.g., 5J]). The graphs G and G'^'sn j-Q^y include all states of the system or a 
fraction of them with the rest summed over, hence the level of the coarse graining can be 
chosen. For example, if one is interested only on the molecule internal dynamics, summing 
over all states of the potential may turn out useful. For a RD-type model this would mean 
that the dimension of the graph is reduced by a factor of 1/SL, which also makes the 
numerical optimization easier. 

The dominating processes simply give a collection of the most probable transitions that 
the molecule can go through successively, thus giving information about the types of pro- 
cesses that are important. The probability for the (complex) molecule to precisely follow 
such fixed paths is of course very small. Because of this, it would be hard and time consum- 
ing to try to identify dominating processes from the simulation or experimental data. Our 
proposed graph analysis is simple and can in principle be done for all finite discrete stochas- 
tic non-equilibrium systems which have non-zero currents. Whether this analysis is worth 
the effort (i.e., if 7 does contain interesting information), depends on the complexity of the 
system and the importance of the molecule internal dynamics to the transport process. 



B. Motor efficiency 

The efficiency of the molecular motor is an important aspect, especially for non-artificial 
molecular motors that have limited energy available. In the literature, there are several 
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definitions of tlie efficiency for Brownian motors, see, e.g., Refs. jll. Issl 59| . Here we adopt 
tlie basic t tier mo dynamic definition that relies on tlie constant load force F on the polymer, 
which means that the output power of the motor is vF. The input power H^in comes from 
turning the potential on, thus forcing the polymer periodically in a higher energy state 
depending on its location. This approach is different from the model where the molecule 
gains constant amount of energy by, e.g., ATP hydrolysis. We assume that the energy is 
dissipated, when the polymer goes back to lower energy state, i.e., this energy is not taken 
into account by reducing it from the input energy. By assuming that transitions between 
potentials of type 1 system are cyclic (i.e., Vi V2 —)■... V5 —)■ Vi), the input power for 
stochastic and deterministic potential schemes can be written as 



Ef=i Ee. max [0, E,+i(e) - E,{e)] T;^P{es) , type 1 



dE{t,e) 
' dt 



P{e,t) 



types 2-3 



where e and every include L3^~^ states. Since the type 2 potential has discontinuities 
in t, one can define dE{t,e)/dt := Ef=i [-^s+i(e) — -^s(e)] S{t — Efc=i^fc)- The efficiency is 
defined by 77 = 

Although the efficiency of the flashing ratchet model is very low for single particles (see. 



recent work 



e.g., Ref. 56[), it can be greatly increased for some many particle systems as shown in the 
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13|. Besides the efficiency, we are also interested in the stopping force F^top 



which, when applied, causes the average drift go to zero. It is expected that the stopping 
force gets larger as increases, as seen in Ref. js]. 



C. Algorithms 

When dealing with large linear systems (of the order of 10^ states and beyond), one must 
really pay attention to the convergence properties and therefore the choice of the numerical 
methods are important. In this Paper we have three types of numerical problems to solve 
-Pstat- For the fully stochastic system (type 1) we used the Arnoldi and bi-gradient stabi- 
lized (BiGradStab) methods (drift and diffusion), for on/off deterministic system (type 2) 
adaptive Runge-Kutta 4-5 method and for smooth continuously deterministic system (type 
3) quasi-minimal residual (QMR) method. The solution of the type 1 problem is a straight- 
forward eigenstate computation, the other two are more involved integration problems. All 
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computations were performed in Matlab with a modern desktop computer. Solving sta- 
tionary states for the type 2 and 3 potentials were the most time consuming parts of the 
computations. 

When solving the stationary state for type 1, a random initial vector is good enough 
choice, but for types 2 and 3 this is not the case. A better initial guess is needed to 
reduce the computation time. We found that the stationary state of the mean-field operator 
{Hmf — X^feS^feiffe with Xk — Tk/T) is easy to compute and a good one to begin with. In 
many cases, previous solutions can be also used (e.g., when varying T). A random initial 
state however serves as a good check of the numerics, since the results must not depend on 
the choice of the initial state. 

The stationary solution for the type 3 can be found with the same manner as for the 
type 2 (RK45), which however requires that the operator is available for all t e [0,T] and 
are either re-build every step or loaded from the memory. The other way (which we used) 
is to solve the larger linear equation problem as a first order discretization in time, 

XX./ X P(t + At)-P(t-At) 

mm «^ ^ — — (4) 

where At = T/M, M being the number of discretization steps. We found that M — 30. ..60 
is accurate enough. In the matrix form this leads to the problem HP — A, where H includes 
H{t) for all M time-steps and the discretization operator, and the normalization is preserved 
with Ai — Vii — LY, 2LY, MLY otherwise zero. As before, the time-dependent diffusion 
coefficient is found by solving another hnear problem. For these linear systems the QMR 
method turned out to be well converging (LSQR is also a fool-proof method, but very slow). 

To maximize or minimize the velocity ^(q, T) for charges and the temporal period, non- 
linear optimization can be carried out with the standard sequential quadratic programming 
method. To find the velocity, the generator if(q, T) must be constructed several hun- 
dred/thousand times because of changes in the transition rates. Efficient implementation 
presumes that this process is fast, which is achieved for example by manipulating the re- 
quired matrix elements directly in the memory instead of re-building the whole matrix. The 
choice of the initial state is crucial (as usual for optimization problems) and a random state 
is used with several repetitions to confirm the global optimal point. A symmetric initial 
charge distribution easily leads to a local optimal point with a symmetric charge distribu- 
tion (as seen in Section IV C). If q is fixed, optimization can be replaced by interpolation, 

12 



since function v{T) is very smooth. 

The best known exact algorithms to find the optimal cycle ratio have the complexity 
0{nm) 60| , where n and m are number of vertices and edges, but in practice these algorithms 



are not the fastest ones 



541 • We applied an improved version of the Howard's method 61 1 



implemented in the Boost C++ library. There also exist brute-force methods to efficiently 
find (enumerate) all cycles in graphs 62|, but this approach is limited to very small networks 
and/or cycle lengths. We also tested a simple greedy algorithm where we begin from a single 
edge with the largest weight and start to grow the path by always choosing the edge with the 
largest weight available at the moment, until the path form a cycle (i.e., crosses itself). This 
method however works poorly and an optimal solution is found only for very simple cases 
(e.g., a polymer in strong static field), where the results are also easy to guess beforehand. In 
general situations, the optimal path contains transitions that cannot be chosen by a simple 
greedy algorithm. 



IV. RESULTS FOR THE DIFFERENT POTENTIAL AND POLYMER TYPES 

Since both RD and FM models include a large number of parameters, some of them must 
be fixed, primarily those that have a minimal qualitative impact on the results. In addition 
to (reptons), other parameters in the models have the following interpretations: 

• The environment -H- the potential V{x, t) = V{x + L,t + T) 

• The medium <(-)■ tube deformation Q (0 for RD, 1 for FM) 

• The polymer internal fine-structure -f-?- charges in q 

The single most important parameter is the period T of the potential, which is also one of 
the easiest one to control in experimental set-ups. The parameter Q models the porosity and 
viscosity of the medium by either restricting polymer strictly into the reptation tube {Q = 0) 
or not (fi = 1). As before in Ref. j28|, we set S* = 2 and L = 3 to achieve a both maximal 
N/L ratio and keep feasible matrix sizes. The fiashing ratchet is Vi(l) = Knax, Vi(2) = 
Knax/2, Vi(3) = and V2(x) = Vx, and the traveling potential Vi(l) = ¥2(2) = Vmax and 



there is an illustration of these 



zero for 1/i(2), 1/i(3), 1/2(1) and 1/2(3). In Fig. 2 of Ref. 
potentials. Time symmetry parameter x = Ti/T is fixed to 1/2 for the fiashing ratchet 



13 



potential and 1/4 for the traveling one. The maximum potential strength Knax has only 



a small effect on the results and is set to unity (with one exception in Fig. [7]) 66|]. The 
direction of the potentials is set up in such way that the expected "main drift" is always 
positive and the inverse drift (if present) is negative. 

With the definitions in Section II, we study the following three types of time dependent 
potentials 

• Type 1: stochastic on/off switching 

• Type 2: deterministic on/off switching 

• Type 3: deterministic smooth cosine-type modulation. 



A. Comparison of time-dependency schemes 

First we compare the differences of the potential time-dependency schemes in the flashing 
ratchet potential, for which the differences are more distinct. In Fig. [3] we have plotted v 
and Pe of = 5 and 9 (similar behavior is observed for other values of A^) RD and FM 
polymers as a function of T for all three time-dependency schemes. 

Some clear differences between the schemes can be seen. The maxima for the drift and 
the Peclet numbers are reached for smaller T for type 1 than for types 2 and 3. The type 
2 scheme has the largest v and type 3 the smallest, and the same goes for Pe. However, 
this order changes for the inverse drifts, where types 2 and 3 are equally good. The time- 
dependency scheme turns out to have an effect on the current inversion phenomena, since 
the type 3 scheme is able to invert all RD polymers with N > 2, whereas types 1 and 2 
only those with > 5. Despite this, the differences between types 2 and 3 are small (type 
2 being slightly "better") and we now concentrate only on types 1 and 2. 



B. Motor efficiency and stopping force 

In Figs, m and ID we show the maximum efficiency r/max = max r7(T) of the RD and FM 
polymers as a function of a load force F = YliLi ^^i^ where E is the field strength, with 
flashing ratchets and traveling potentials of the type 1 and 2. The points where ?7max(-^) = 
for F > define the stopping forces Fstop- Insets of the figures show the same data scaled 

14 




FIG. 3: (Color online) Drift velocity and Peclet number for 5 (blue (gray) lines) and 9-repton 
(black lines) RD (left) and FM polymers (right), for type 2 (solid lines), type 3 (dashed lines), and 
type 1 (dash-dotted lines) schemes. 

with F* = F/Fstop and r]^^^ = ''7max(-^*)/ niax?7jnax(-^*) for each polymer size, which reveal 
the shapes of the curves. 

We notice that for the FM polymers the efficiency is generally larger and they can main- 
tain their drift in an opposing field better than the RD polymers in the ratchet. When 
plotted as a function of E, there is a constant stopping field for all > 3 FM polymers in 
both potentials with values around -0.0026/-0.0016 for type 1 and -0.0038/-0.0043 for 
type 2 ratchet/traveling potentials. This results from the fact that the reptons of the FM 
polymer are less correlated than those of the RD polymer and the FM polymer thus behave 
more independently . For the ratchet, the type 2 scheme is found to be 2-4 times more 
efficient and can withstand almost double load force when compared with the type 1. The 
stopping force is larger for FM polymers. For the traveling potential, differences are more 
drastic, as for the type 2 scheme the stopping force is about two times and the efficiency 
almost one order of magnitude larger when compared to the type 1 scheme. Rescaled curves 
reveal that despite the large differences in scales, shape of the curves are almost identical 
for all polymer lengths and both types. 

The numerical values of the efficiency are very small. This is a generally known trait 
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FIG. 4: (Color online) Maximum efficiency for RD (black solid lines) and FM (blue dashed lines) 
polymers as a function of the load force F with = 1...9 in type 1 (left) and type 2 (right) flashing 
ratchets. In each case, the rightmost curve is for = 9 and the bold lines (the less interesting 
special cases = 1, 2) are shared for both RD and FM polymers. Insets: Rescaled data ^max ^ ^ 
function of F*, with black triangles for RD polymers and blue squares for FM polymers. 




FIG. 5: (Color online) Maximum efficiency for RD (black solid) and FM (blue dashed) polymers as 
a function of the load force F with = 2...9 in type 1 (left) and type 2 (right) traveling potentials. 
In each case, the rightmost curve is for A^ = 9 and the bold lines (the less interesting special case 
A^ = 2) are shared for both RD and FM polymers. Insets: Rescaled data r/j^ax a function of F* , 
with black triangles for RD polymers and blue squares for FM polymers. 

especially for flashing ratchet models 56|, but it also results from the choice of the rates, 
since the velocity plays dominating role for the efficiency. By the use of the optimized 
parameters (e.g., Knaxja^^q), efficiency could be increased by couple orders of magnitude. 
Results show that Fstop increases as a function of A^, which is in agreement with some 
previous work B Q. TKe .owe.. ..... a. po,.. .e. .o... fo. 
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all other but the type 1 traveling ratchet, which is surprising. 
C. Non- uniform charge distributions 

Extensive computations were carried out to find the charge configurations with the largest 
possible V in forward and backward transport and Pe for various polymers and parameters. 
It was found that changes in the drift are so large that one can safely limit to maximizing 
V alone, since in this case Pe is dominated by the drift. In the following, some of the 
optimization results are presented for the 8-repton polymers in the type 1 potentials. The 
basic model with an uniform charge distribution (g^ = IVz) is also shown for comparison. 

In Fig. [6l the properties of the RD polymer in the fiashing ratchet are plotted as a func- 
tion of T with configurations that give maximum drifts for positive (forward) and negative 
(backward) directions, and the corresponding optimization results are called either positive 
or negative. We found that the positive direction is always maximized by putting all charge 
near either of the heads, but charging the head reptons does not necessarily lead to the 
largest current. This holds for both RD and FM polymers for all studied polymer lengths 
up to = 13 at least. In this situation only one repton feels the potential and very large 
transition rates are generated by the exponential function (see Sec. II B). This one repton 
then forces the whole polymer to advance. 

The optimal charge distributions in the negative direction are more interesting, since the 
large accumulations of the charge are not seen and the charge is distributed over several 
reptons. Symmetric distributions results that neither of the heads are leading and are 
forced to compete with each other. This would be very inefficient in constant-field transport. 
Repeating the optimization computation several times, additional distributions very close to 
the first one are found. Similar local optima are also found in other cases, which complicate 
the search for the global optimal distribution. This is demonstrated in Fig. [71 where we fix 
T = exp(5.5) and show the drifts given by the three local optimal distributions as a function 
Vmax for the RD polymer in the fiashing ratchet. At Knax ~ 1.126 the non-symmetric 
distribution becomes the fastest one. 

In Figs. [8] and [9] we show the same analysis for the traveling potential. For RD poly- 
mers, the optimal distributions have no "clear" or symmetric structure, only near optimum 
symmetric distributions are found. As seen above for the fiashing ratchet, the negatively 
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FIG. 6: (Color online) Drift velocity, diffusion coefficient and Peclet number for the 8-repton RD 
polymer in the flashing ratchet with uniform (blue dash-dotted lines), negatively optimized (red 
dashed lines), and positively optimized (black solid lines) charge distributions as a function of the 
temporal period T. The histogram shows the charge distribution along the polymer for each case 
in the same order. In the leftmost figures v and Pe for the positively optimized case have been 
scaled by an additional factor 1/10. 

optimized polymers are actually faster in both directions when compared with the uniformly 
charged polymers. 

In conclusion, the charge distribution has a large effect on the polymer transport velocity 
and coherence on the flashing and traveling potentials. Since the drifts generated by the 
ratchet effect are generally very small and difficult to observe, this could be of interest from 
the point of view of applications. In the next Section we show that different distributions 
also lead to different kinds of transport mechanisms. 
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FIG. 7: (Color online) Left panel: Drifts generated by three locally optimal charge distributions 
of the 8-repton RD polymer to negative current direction as a function of l^nax- Right panel: The 
charge distributions at Vmax = 1-1734, with the leftmost bar corresponding to the dash-dotted blue 
line, the bar in the middle to the dashed red line, and the rightmost bar to the solid black line. 




FIG. 8: (Color online) Drift velocity, diffusion coefficient and Peclet number for the 8-repton RD 
polymer in the traveling potential with uniform (blue dash-dotted lines), negatively optimized (red 
dashed lines), and positively optimized (black solid lines) charge distributions. The histogram 
shows the charge distribution for each case in the same order. 
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FIG. 9: (Color online) Drift velocity, diffusion coefficient and Peclet number for the 8-repton FM 
polymer in the traveling potential with uniform (blue dash-dotted lines), negatively optimized (red 
dashed lines), and positively optimized (black solid lines) charge distributions. The histogram 
shows the charge distribution for each case in the same order. 

V. RESULTS FOR THE INTERNAL DYNAMICS OF THE POLYMERS 
A. Time evolution of the observables 

To gain better insight in the internal dynamics of the polymer we now turn to the expected 
values of the four observables Z (zero-bond), G (total length), K (kinks) and H (head-to- 
head distance) for the RD polymer. In Figs. [10] and [H] we have plotted the stationary 
state time-evolution of the observables against each other with the 8-repton RD polymer in 
the type 2 ratchet and traveling potentials with several values of T. The previously found 
optimized charge distributions are used. Note that these distributions are only approximately 
optimal for the type 2 potentials, but this approximation is found to be very good. For small 
T, the observables are near their mean-field values (large spots in the figures), which are 
independent of t. For very large T, the curves "freeze" (bold lines) since the stationary 
states are reached before the potential is switched. 
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FIG. 10: (Color online) Time evolution of the 8-repton RD polymer observables in the type 2 
flashing ratchet with uniform (red (gray) lines), positively optimized (black lines) and negatively 
optimized (blue (light gray) lines) charge distributions. In each case, the big spots correspond to 
the (mean-field) limit T — > 0, the bold curves show the T — >• oo limit and the other curves the 
behavior for a few selected finite values of T. 

For the ratchet, the maximum positive current (black lines) is a result of small changes in 
the polymer average shape, which is caused by the fact that only a single near-head repton 
is charged and the rest of the polymer is in pure random motion. The maximum negative 
current (blue (light gray) lines) however is a result of more complex processes, which cause 
much more variation in the average shape, even more than for a polymer with uniformly 
distributed charges (red (gray) lines) with all the reptons charged. There is almost one-to- 
one correlation between G — H pair (as expected), which results that the phase trajectories 
for the G — Z and H — Z pairs are almost indistinguishable, and therefore the pairs G — H 
and G — Z are not presented here. The connections between other pairs are more involved, 
especially between H—K and G—K. For them, the current inversion is seen as a deformation 
of loops between K — Z, H — K and G — K pairs for uniform and negatively optimized 
polymers (no current inversion for positively optimized polymer). 

For the traveling potential, the curves are more distinguished from each other and are 



21 




FIG. 11: (Color online) Time evolution of the 8-repton RD polymer observables in the type 2 
traveling potential with uniform (red (gray) lines), positively optimized (black lines) and negatively 
optimized (blue (light gray) lines) charge distributions. In each case, the big spots correspond to 
the (mean- field) limit T — > 0, the bold curves show the T — >• oo limit and the other curves the 
behavior for a few selected finite values of T. 

more complicated. There are fast deformations in the curves as the time goes on. There is 
a clear similarity between Figs. [10] and [TT] Positively optimized polymers have the smallest 
spread in the observables and negatively optimized the largest. This is similar behavior 
as seen for the flashing ratchet, albeit the potential and the charge distributions are very 
different. The results show that there is a clear connection between the average polymer 
drift magnitude and direction, and shape deformations. Deformations during ratcheting for 
one's part depend strongly on the charge distributions. 

In Figs. [12] and US] we have plotted the relaxation of the observables in the flashing ratchet 
and the traveling potential for the 8-repton RD polymer with uniform and optimized charge 
distributions. The data is the same as shown in Figs. [10] and [11] for the large T limit (bold 
lines). For the observables, the largest changes are observed in roughly the same time scale, 
around ln(t) ^ 2. Stationary values for observables for the positively optimized polymer in 
the flashing ratchet are independent of the potential state (on or off). In addition to the 
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FIG. 12: (Color online) Relaxation in real time t of the 8-repton RD polymer in the type 2 ratchet 
potential, with uniform (red (gray) lines), positively (black lines) and negatively (blue (light gray) 
lines) optimized charge distributions. Dashed lines (when present) are for the "on — )> off" and solid 
ones for the "off — )• on" processes. 

kink dynamics, large differences are shown in zero-bond dynamics. Note that for positively 
optimized polymer, values remain unchanged during "on — )■ off" switching and are therefore 
not shown in the figure. This is because, in the steady state, the potential has no effect on 
the conformations of the polymer, which would require more than one charged reptons. For 
the traveling potential, the time-evolution of the observables is more complex. 

B. Network analysis 

To further understand the formation of the net drift, we now turn to the network analysis 
of the steady state currents. We concentrate on the RD polymer of the type 1 in the flashing 
ratchet and the traveling potential with uniform and optimized charge distributions. The 
temporal periods T are chosen such that they result in the maximum current (4 values of T 
for both potential types). The graphs C^'S"^ containing the steady-state net currents between 
the states are then computed. We have summed over all the potential states {SL degrees of 
freedom) so that only the internal states of the polymer remain. After these steps we have 
eight different graphs with 5832 non-zero directed edges in each of them. 

Let us flrst analyze these G^^^^ graphs by defining the arrays S with elements Si [i = 
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FIG. 13: (Color online) Relaxation in real time t of the 8-repton RD polymer in the type 2 traveling 
potential, with uniform (red (gray) lines), positively (black lines) and negatively (blue (light gray) 
lines) optimized charge distributions. 

1...5832) that contain all edge weights of the graphs G^^^^ in an increasing order. In Fig. [13] 
we show Si and their cumulative sums. The total drift v is then recovered as the sum of 
all Si (c.f. Figs. 6 and 8) and a non-zero drift is produced when the S curves are tilted 
due to the external forces. Although the system is far from equilibrium, only a slight tilt is 
observed and there are no single dominating edges. 

We now turn to the dominating transport cycles of the polymer motion by analyzing 
the paths in G^^s". This results in cycles with lengths of the order of 10. It is found that 
the common transportation type is such that we call "si-S2-scheme" consisting of cyclically 
accumulated (lengths Si and S2 with |si — S2I = 1) and elongated parts of the polymer. 
Corresponding to the direction of moves, this scheme can be either positive (up) or negative 
(down). To illustrate the scheme, we have sketched the positive 4-5 scheme in Fig. [151 The 
numbered arrows indicate the order and direction of the corresponding repton moves. After 
all marked moves are done, the initial state is recovered and the cycle is repeated. In the 
five situations out of eight studied here, the dominating cycle is the si-S2-scheme. 

In Fig. [T6]we show the remaining three situations that are not of the type above. Note that 
for negative transport in the ratchet with the uniform charge distribution, the mechanism 
is almost the negative 4-5-scheme. 

In Table [Tj we have collected the core results of this Section. For comparison, there 



24 




FIG. 14: (Color online) Upper panels: Ordered elements Si of the matrix G^^^^ for the 8-repton 
RD polymer in the type 1 flashing ratchet (left column) and traveling potential (right column) 
to positive (black lines) and negative directions (red (gray) lines), with uniform charge distribu- 
tions (solid lines) and optimized distributions (dashed lines). Lower panels: The corresponding 
cumulative sums Yl\=i ^k- 




FIG. 15: Illustration of the positive 4-5-scheme for the 8-repton polymer. The arrows and numbers 
indicate the direction and the order of the transition for the corresponding reptons. For clarification, 
the process is shown here in two parts. 

is a ratio of the average cycle drift Vc, divided by the average drift by the all transitions 
^aii = f/5832 in the last column. This ratio is significantly larger for uniformly charged 
polymers, indicating that the optimization process increases the drifts for large number of 
paths and makes differences between paths smaller. It is also somewhat surprising that there 
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FIG. 16: Dominating cycles for the backward drift in the travehng potential with the uniform charge 
distribution (a), the forward drift in the traveling potential with the uniform charge distribution 
(b), and the backward drift in the ratchet with the optimal charge distribution (c). 

is not much difference between the leading mechanisms for forward or backward motion, and 
for uniformly charged polymer in ratchet it is actually the same. One can therefore conclude 
that the current inversion for the RD model is not caused by some abrupt 'phase transition', 
but gradual changes in the probability distribution along internal states. 

We carried out a similar analysis also for the full system without summing over S and 
L, in which case cycles have up to 30 states and there are some modifications to the pure 
Si — S2 schemes. However, these cycles are too lengthy to be reported here. It was found 
that sometimes summing over the potential states is necessary to find a non-stationary cycle 
and sometimes the summing leads to stationary cycle. 

VI. DISCUSSION 

We have analyzed the properties of Rubinstein-Duke polymers with some modifications, 
including tube breaking and non-uniform charge distributions, in time-dependent potentials. 
The aim of this work was to further study the properties of complex molecules in out-of- 
equilibrium conditions and especially the ratchet effect. 
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Case Cycle 


Vcl V^ll 


Ratchet potential 


unif. F pos. 4-5 


46,4 


unif. B neg. 4-5 


554 


optim. F pos. 2-3 


6,5 


optim. B Fig. [T6lfc) 


164 


Traveling potential 



unif. F Fig. HSl^b) 60,0 

unif. B Fig. HS^a) 371 

optim. F pos. 2-3 29,4 

optim. B neg. 1-2 91,3 

TABLE I: Dominating cycle types for polymers in ratchet and traveling potentials for forward (F) 
and backward (B) transport, with uniform (unif.) and optimized (optim.) charge distributions. 



In the first part of the study, we extended the previous work reported in Ref. [28[ by con- 
sidering deterministic ratcheting mechanisms, the energetic efficiency and optimized charge 
distributions of the polymers. It was found that the deterministically flashing potential is 
superior when compared to a smoothly varying and stochastic potential for velocity, co- 
herence and efficiency. However, despite "scaling differences" in drift and diffusion, the 
time-dependency scheme seems to have a minor effect on the qualitative results. By us- 
ing the stochastic scheme, we computed the optimal charge distributions to maximize the 
steady-state velocity in flashing ratchets and traveling potentials. The differences between 
these and the uniformly charged polymers were found to be drastic. Changing the charge 
distribution also changes the mechanism of how the polymer reshapes itself with respect to 
the potential. 

In the second part, the current inversion phenomenon was investigated in detail by us- 
ing the optimal charge distributions. The expected values of certain macroscopic observ- 
ables (e.g., length and zero-bond count) were computed and large differences between dif- 
ferently charged polymers were found. To find how the polymer actually moves in the 
non-equilibrium steady state, we proposed a simple graph analysis method to find most 
probable series of state transitions (=path) based on the probability currents. For a pe- 
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riodic system such a path is found as a solution of the optimal cycle ratio problem. This 
method is suitable in situations where a huge network is generated by some automated fash- 
ion or measurements and cannot be analyzed "manually" (e.g., Kinesin network in Ref. {gS]). 
This method was then used to identify the dominating processes of the polymer transport 
and was found to be very useful to piece together polymer motion. However, the general 
usefulness of this analysis depends on the model and it would be of interest to test it for 
other complex out-of-equilibrium systems and also with non-periodic boundary conditions. 
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Appendix A: Details of the operators 

In this Appendix the polymer state operators are explained in more detail with some 
practical aspects of constructing them. 

1. Definition of Eq. (1) 



The explicit definitions of the operators in Eq. (1) are as in Ref. 28|], but due to the 
arbitrary charge distributions there is an additional charge dependency in the functions L 
and R: 

Ai{q) = {Riq, I) + L(/)}n0,i,i - - L{q, l)^^ 

+ L(g, l)hA,i,i - L{q, l)di^i + R{q, l)nBxi - R{<1, OKi 

By^iiq) = {R{q, l + f{N- 1, y)) + L{q, l + f{N- 1, y))}n^,N-i,y,i 

- Riq, l + f{N- 1, y))aj,_,_^ ^ - / + /(iV - 1, l/))&jv_i,,,, 
+ L{q, I + f{N - 1, y))nA,N-i,y,i - L{q, I + /{N - l,y))aN_i^y^i 

+ R{q, l + f{N- 1, y))nB,N-i,y,i - R{q, l + f{N- 1, y))bN-i,y,i 
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Mi^y,i{q) = R{q, I + f{h y)){nA,i,y,in0,i+i,y,i + nis,i,y,inB,i+i,y,i - ai,y,ial+i^y^i - b\y,fii+i,y,i) 
+ L{q, I + /(i, y)){n0^i^y^inA,i+i,y,i + ns^i^y^in^^ij^i^y^i — 0'l^y^iO'i+i,y,i — h,y,ib\j^i y j) 

,i,y,l1T'B,i+l,y,l + ^0,i,j/,«^0,t+l,|/,« <^i,y,lbi+l,y,l bi,y,fii+l,y,l) 

+ VLL{q, I + /(i, y)){nB,i,y,inA,i+i,y,i + - 6i,y,zai+i,y,z - 4„/J+i „ ,), 



where f2 = for RD polymers and f2 = 1 for FM polymers, and 

nz,i,i = nin^^i 
/i-i 

^i,y,i ~ I W ^g(yd)d 
^z,i,y,i = ni ^a{y,j),j 



with X G {a, 6, a^, 6^}, z G {A, 0, B}. The function 5f(y, i) G {A, 0, 5} gives the state of the 
zth bond in the configuration and the function / 

i 

f{i, y) = ^{^y\nA,i - nB,i\^y), l<i<N-l 
i=i 

gives the position of the repton i + 1 in marker-centered coordinates. The detailed forms 
of the functions g and / depend on the selection of the state basis. With these functions, 
the formal definitions for the macrostate observables, i.e., zero-bonds, kinks, head-to- head 
length and the total length, of the repton polymer are 

yeFl : 4^{l<z<N;g{y,t) = 0} = n,, 

y e K^^ : < t < N - l;g{y,t) = A/B A g{y, z + = B/A} = nk, 

|/(iV-l,y)|=n,, 
y G : max [/{k, y) - f\l, y)] = n„ A;, / = 1, 2, iV - 1 

where n^, n/i, = 0, 1, iV — 1 and Uk = 0, 1, ...,N — 2. One can verify that i^Ff^ > jj^F^^ 
holds for all i. By using above sets F and equation ([2]), measure operands can be constructed 
and expected values computed. The practical procedure to form all the required operators, 
especially the previous observables, is explained below. 
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2. Operator construction 



Since the stochastic generator and measurement operands used in this work are shghtly 
more complex than in the previous works regarding the RD model, we show in some details 
how the idea of the recursive operator construction work in the current case. Whereas 
small operands can always be build directly, recursive construction is practically a must for 
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64j |. For simplicity. 



large systems and nowadays widely used in DMRG computations 
we concentrate only on (discrete) state measure operators, which in the natural basis are 
diagonal matrices. 

Let |0]^, O*. } be a set of macrostate operators for the system with i sites, which 
includes all the necessary operators that are required when adding a new site. Here site is 
a general term, which for example could mean single particle states for classical systems 
and spin states for quantum systems. By using the usual product state formalism, assume 
that the new sites are added on the right such that |new state) = |old state) ® |new site). 
The basic algorithm to add new sites (until A^) goes as follows 



1. Build an initial set of operand(s) Oy, where y = I, ■■■,yi- 

2. For all m = 2, 3, and y = 1, ym, build: 

{k,j)=K{ym-i,y) 

3. Build the full operand = YJijli ^yO^ , 

where ym is the total number of operands required for the size m system. The details 
of how to construct a new set of state operators for the enlarged system by joining the 
states of the new site and the old operators are hidden in the function K{ym-i,y)- The 
complexity of this function and the number of required operators ym depends on the type 
of the operand. Practically it is the ym that determines the computational effort needed to 
build large operators, since K is more or less just keeping book of how to join operators. 

We now concentrate on the RD-type model for which sites mean bond states between 
the reptons. The polymer state operators have the following values for ym, given with brief 
explanations 
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• Zero-bonds: |/m = m + 1 (number of zero-bonds) 

• Kinks: Um = max {3, 3 (m — 1)} (number of kinks and state of the rightmost bond) 

• Head-to-head length: = 2m + 1 (signed distance between the heads) 

• Total length: = (1 - /)(3/ - 2m - 7) > (1 + m) (3 + m) /3, where I = \^] (see 
example below) 

The number of required operands is therefore oc m^ for total length and oc m for others. 

We now consider a concrete example for a total length operator, which is the most complex 
operator used in this paper. When one enlarges the size of this operator with new particles, 
one must keep track of the maximum distances of the rightmost repton from all the other 
reptons. For example, in Fig. [1] these distances would be 2 (from repton nr. 4) and (no 
any reptons below). We define these as up (u) and down (d) distances. Total distance is 
then d + u. 

In Fig. [T7] we show all 9 microstates of the 3-repton polymer. Since there are five com- 
binations for u and d distances, the macrostate operators o)' 2)' i) formed 
with each of them including one or more microstates. This is shown in the figure with 
red numbers in {d, M)-plane. When a new repton is added, function K is used to combine 
old macrostate operators with state operators of the new site ({n^, ""-0, n^}) and hence the 
number of macrostate operators is increased by three. Examples of the required operations 
includes O^^ 0) ~ ^fi o) ® ^{1 1) ~ 0) ® ^-B- After addition, there are 27 microstates 

in eight macrostate operators (blue numbers in the figure). Note that in the actual com- 
putations only the information about the d and u values is needed. Here the tracking of 
the microstates was done for illustration purposes only. As more reptons are added, the 
"triangle" that presents available {d, u) states gets larger. 
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FIG. 17: (Color online) All 9 configurations (microstates) of the 3-repton polymer with 5 macrostate 
operators that are formed from them. The macrostates are indexed by the corresponding (d, u). As 
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This is clear, since in the linear response regime only the total force directed to the object 
counts and the diffusion coefficient is given by the Nernst-Einstein relation. Similar equivalence 
does not hold for systems involving non-homogeneous potentials. 

Note that the choice V^ax = 1 in this paper is equivalent with Vmax = 1/2 of the previous 
paper, where the factor 1/2 was dropped in the definition of rates Hij 
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